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BACKGROUND OF THE INVENTION 

1 . Field of the Invention 

The invention relates generally to the field of imaging. More particularly, the 
invention relates to methods and systems for image reconstruction using cameras, such as 
Compton cameras. 

2. Discussion of the Related Art 

In vivo imaging methods, such as Single Positron Emitting Computed Tomography 
(SPECT), use gamma radio tracers to track biochemical, molecular, and/or physiopathological 
processes of various human diseases. Further, such imaging methods may provide detection 
of contaminates in a nuclear facility and nuclear waste sites as well as serve as a defense 
utility by imaging radiation activity on missiles, planes, etc. Radiation emitted from the target 
of interest are detected by a gamma ray camera device of the imaging system which forms an 
image of the target based on the concentration and distribution of the radioactive material, 
e.g., gamma radio tracers within the target. 

Conventional gamma ray cameras usually include a plurality of detectors For 
example, the Anger camera includes collimators within the detectors to limit radiation 
trajectories observed by the detector. In addition, conventional gamma cameras are stationary 
within an imaging system, allowing only one view of the target to be observed. However, 
conventional cameras are inefficient because the cameras lack the ability to image two or 
more isotopes at the same time. 

Recently, Compton cameras have been integrated into imaging systems as an 
alternative to the conventional gamma cameras to further enhance and improve upon the 
quality of images being reconstructed. For example, the Compton camera can detect one to 
two orders more radiation emissions, such as photons emitted from a target, than a 
conventional camera. In addition, the Compton camera can readily image a relatively wide 
range of energies. Generally, Compton cameras include two semiconductor detectors 
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configured in parallel with one another. The first detector may be capable of measuring the 
radiation emitted from a target, such as a photon emission. In particular, the first detector 
measures the point in which the photon contacts the first detector as well as the amount of 
energy lost by the photon when the photon goes through Compton scattering within the 
detector. As a result of the scatter, the photon travels in a new direction and interacts with the 
second detector, in which the second detector can measure the point at which the photon 
contacts the second detector. 

Current methods for image reconstruction utilizes the cone projection data collected 
from the Compton camera and convert the data into plane projection data. Such methods may 
include using an infinite series expansion to calculate the plane projection data. However, 
such methods require an a large number of calculations to be performed, and therefore 
requires a greater number of resources to be utilized and may produce poor quality images. 

The referenced shortcomings are not intended to be exhaustive, but rather are among 
many that tend to impair the effectiveness of previously known techniques concerning image 
reconstruction however, those mentioned here are sufficient to demonstrate that methodology 
appearing in the art have not been altogether satisfactory and that a significant need exists for 
the techniques described and claimed in this disclosure. 

SUMMARY OF THE INVENTION 

Thus, there is a need for methods and systems that provide an efficient and accurate 
techniques for image reconstruction. 

In one embodiment, the invention involves a method. The trajectory of a photon from 
an object through a first detector determines an apex of the cone at the point of intersection 
with the first detector. The trajectory of the photon through the first detector onto a second 
detector determines an axis of symmetry of the cone. Using the apex and axis of symmetry of 
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the cone, a set of integrals are calculated and used for image reconstruction. 

In another respect, a method for image reconstruction. A set of conical integrals, such 
as surface integrals or integrated cone-beam line integrals are calculated to satisfy a 
completeness condition. The set of integrals are then related to a distribution of radioactivity. 

In another aspect, a computer readable medium including instructions is used to 
determine an apex and an axis of symmetry of a cone. The apex and the axis of symmetry is 
use to calculate a set of integrals, such as surface or integral line integrals, to satisfy a 
complete condition. After the calculation of the set of integrals are complete, instructions are 
provided to reconstruct an image. 

In another respect, a system may include a camera, such as a Compton camera and 
first and second detectors that are configured to obtain data to satisfy the completeness 
condition. 



These, and other, embodiments of the invention will be better appreciated and 
understood when considered in conjunction with the following description and the 
accompanying drawings. It should be understood, however, that the following description, 
while indicating various embodiments of the invention and numerous specific details thereof, 
is given by way of illustration and not of limitation. Many substitutions, modifications, 
additions and/or rearrangements may be made within the scope of the invention without 
departing from the spirit thereof, and the invention includes all such substitutions, 
modifications, additions and/or rearrangements. 

BRIEF DESCRIPTION OF THE DRAWINGS 
The drawings accompanying and forming part of this specification are included to 
depict certain aspects of the invention. A clearer conception of the invention, and of the 
components and operation of systems provided with the invention, will become more readily 
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apparent by referring to the exemplary, and therefore nonlimiting, embodiments illustrated in 
the drawings, wherein like reference numerals (if they occur in more than one view) designate 
the same or similar detectors. The invention may be better understood by reference to one or 
more of these drawings in combination with the description presented herein. It should be 
noted that the features illustrated in the drawings are not necessarily drawn to scale. 

FIG. 1 illustrates an image reconstruction system, according to embodiments of the 
disclosure. 

FIG. 2 is a block diagram of components within a Compton camera, according to 
embodiments of this disclosure. 

FIG. 3 illustrates an implementation of a completeness condition, according to 
embodiments of this disclosure. 

FIG. 4 shows a geometry for defining probability functions of a two-step 
reconstruction method, according to embodiments of this disclosure. 

DESCRIPTION OF PREFERRED EMBODIMENTS 
The invention and the various features and advantageous details thereof are explained 
more fully with reference to the nonlimiting embodiments that are illustrated in the 
accompanying drawings and detailed in the following description. Descriptions of well known 
starting materials, processing techniques, components and equipment are omitted so as not to 
unnecessarily obscure the invention in detail. It should be understood, however, that the 
detailed description and the specific examples, while indicating preferred embodiments of the 
invention, are given by way of illustration only and not by way of limitation. Various 
substitutions, modifications, additions and/or rearrangements within the spirit and/or scope of 
the underlying inventive concept will become apparent to those skilled in the art from this 
disclosure. 
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In one embodiment, discloses a reconstruction method which calculates a finite set of 
data to obtain a relationship between the calculations and a distribution of radioactivity from 
an object is provided. In one embodiment, the object may be a human or an animal having 
radioisotopes delivered internally for imaging. In another embodiment, the object may be a 
nuclear facility or a nuclear waste site in which the reconstruction method may determine the 
distribution of contaminates within the facility. In yet another embodiment, the object may be 
a missile, where a number of nuclear warheads on the missile may be determined. 

To reconstruct an image of an object, such as patient 16, a system 14 may include a 
Compton camera, which receives data from the radioactivity emitting from the object, as 
illustrated in FIG. 1. A controller 12, coupled to the system 14, obtains the data, processes the 
data to obtain desired image(s), and outputs the final image(s) to an output device of choice, 
such as a display monitor 10. In one embodiment, the controller 12 processes a plurality of 
measurements from the camera, such as surface integrals from a cone projection. In another 
embodiment, the controller 12 processes a plurality of integrated line integrals from a cone 
projection. 

I Relating Surface Integrals to a Distribution of Radioactivity 

According to one embodiment of the invention, the image reconstruction may be 
calculated from data collected by a Compton camera. In particular, a sequence of cone surface 
integrals within a predetermined set may be calculated. The predetermined set may identify 
what calculations are needed to reconstruct an image of an object. Subsequently, the sequence 
of surface integrals may be related to the distribution of radioactivity and an image of the 
object may be obtain using Hilbert transforms and the partial derivatives of a three- 
dimensional Radon transform. 

The surface integrals may be calculated from a cone projected from the interaction of a 
photon through detectors of a Compton camera. Referring to FIG. 2, the trajectory of a 
photon 104 from an object may intersect a first detector 100 at point <D. Point O may be 
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defined as an apex of the cone. The photon may undergo Compton scattering causing a 
change trajectory of the photon (as illustrated by arrow 106) through the first detector 100 
onto the second detector 102. In one embodiment, the angle change between the trajectory 
(arrow 104) from the object to the first detector and the trajectory from the first detector 100 
to the second detector 102 (arrow 106) may be defined as a half angle y of the cone where 
vector p, which may extend from point O' through point O, may be axis of symmetry of the 
cone. Thus, symbol S((D, p, vj/) may denote the surface integral of the distribution of 
radioactivity on the one sheet cone whose apex may be the point c£, axis of symmetry may be 
the unit vector p, and half angle may be v|/. 

In order to determine the source of the radioactivity, the distribution of radioactivity at 
point x may be defined as J[ x ): In one embodiment, it may be assumed that f{y) = 0 for 
y > R, where R may be the radius of the distribution activity. The vectors O, p, and % may 
also be described in terms of a global coordinate system. Let vectors <D, p, and % be described 
in terms of a global coordinate system. A local coordinate system, where a "z" axis points in 
the direction of vector p may be used in expressing in a spherical coordinate system by letting: 

A A 

a = a{$, y/) = (cos (f> sin y/ 9 sin </> sin y/, cos y/) r (l) 

where v|/ may be the angle measured from the "z" axis. Using a standard calculus equation the 
surface integral may be calculated, where: 

S(®,/3 9 y/) = siny/ j J/fo + rM 7 a (</>,y/))r dr d<f>, (2) 

where the rotation matrix M T may be defined as 

M r =\p u \/3 12 \/3] (3) 
where p, pn and p l2 are three orthonormal column vectors in <R 3 . 

To relate a surface integral to the distribution of radioactivity, a new function C(j3j) 
may be defined as: 
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Cfo/)=jta \f(p,t)p s (t-t)dt, (4) 

-oo 

where 

„.ritai»i>«. 

[0 otherwise, 

V 

md f{fi>t) m ay be a three-dimensional Radon transform, namely: 

fa ') = J J/fo* + + tPn )dsdt. (5) 

— co— oo 

Accordingly, the function C{/3j) may be the Hilbert transform of the three-dimensional 
Radon transform. (Gel'fand and Shilov, 1964 and Bracewell, 1978). 

The relationship between the surface integrals and the distribution of radioactivity may 
be proven as follows: 

Proof 1 

C{/3, = - lim Js(<D, p, W )p s (cos y) dy. (6) 

0 

Proof: For a function g, 

}f{P,t)g(l)dt= \f(z)g(z-P)dz. (Smith etal., 1977) (7) 

-00 ^3 

By letting g(t) = p £ (o . (3 - i) and making a change to the variables, Equation (7) may be 
written as follows: 

00 

\f{fi,t)p. (*-fi-t)dt= jf(®-z)Ps(x-P)dX- (8) 
Letting % = ra foraeS 2 the right hand side of the above equation becomes: 

+00 

J \f (<D - m) p e (ra • fi)r 2 dr da, (9) 

S 2 * 
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where da may the surface detector on the unit sphere. Evaluating the above spherical integral 
in the same local coordinate system that was used in Equation (2) the above express becomes: 



-KJO/T 2/T 



jj jf(®-rM T a(0,y/)) d<f> p e (r cos y)sin y/ dy/ r 1 dr. 



(10) 



0 0 0 



Now letting z = cos ^results in 

-ho +1 In 

\ \ J/(0 -rM T a(</>, cos"' z))d# p s {rz)dzr 2 dr. 



0-10 



01) 



Taking the limits as e-> 0 of the above and then using Lemma 1 (please see below) results in: 



4-oo/r 2n 



lj jf(®-rM T a(4>,cos- l z)) d</> p £ (z)dz r dr 



(12) 



0 0 0 



Letting * -cos' z and exchanging integrals yields: 

+00 n In 

\ \ ~ r M r a(<f>, y/)) r dr dt/> p E (cos y/) sin y/ dy/. 



000 



(13) 



The conclusion of the Proof 1 may allow the substitution of Equation (13) into Equation (8), 
and taking the limit of Equation (8) as s -> 0 and by making some variable changes. 
Lemma 1 For a function g and a constant k± 0 



Sfe(0/>,(fc) dt = limjjg(t)p £ (t) dt 



(14) 



Proof: Explicitly writing out the left hand side (LHS) of Equation (14) yields: 



LHS = lim 

£-+0 



\g(t)dt 



(15) 



Replacing e / |k| with e results in 
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LHS=lim< 



f-e «, \ 

H 



\gW j 



(16) 



The conclusion of Lemma 1 follows by observing that the right hand side of Equation (16) is 
equal to the left hand side of Equation (14). 

To relate the surface integrals to the distribution of radioactivity a second function 
F(j3j) may also used. For p € S 2 and t € 5R 1 the function is define F{/3j) as: 



F{fi, ifUim 1h £ {£ - iff fat) dt, 



(17) 



where 



p-/or|r|>« 

-1 - II 

■^-/or|f|>* 



(18) 



Za 2 

If the function is known on the set p p ={(j3j): 0 e ~ , \t | < r} then the function of 

fix) can be obtained for all % e 9? 2 (Smith, 1985). 



The relationship between the surface integrals and the distribution of radioactivity is 
completed with the following: 



Proof 2 



(19) 



The proof of Equation 19 follows by performing an integration by parts. (Horn, 1978 
or Smith 1998). 
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As previously mentioned, if the function F(j3j) may be known on the set Pf then the 
function fix) for all e 9? 3 may be found. In theory, if a function is known in an arbitrarily 
small neighborhood of a point then the derivative of the function may be computed at that 
point. As a consequence, the function C{p,i) may be known on the set 

p c ={(fij) 'fi e ~^ >\Z\^R + £f where e is an arbitrary positive number, then the function 

F(pj) may be obtained on the set P F via Proof 2. As such, the function C(pj) may be 

obtained on the set Pc from values S(<D, p, y) of via Proof 1 if the following condition is 
true: 

If for almost every (p,i)sP c there exists an G> such that <D -fi = l where 
S(0, p, v|/) is known for all 0 < y < n. 

For convenience let the symbol P s represent any set of (<D, p, y) that satisfies the above 
condition. Using this symbol, this method of reconstructing the distribution of radioactivity, 
fix) on 5R 3 , from a set of surface integrals, S(0>, p, y), may be representative symbolically 

as: 

5(<D, P,y)onP s Prro/1 ) C (fi, £) on P c *°°' 2 > F(p, t) on P F < 52 ) Smith85 > f( x ) on 5R 3 . 

(20) 

As such, set P s and P c define what integrals are needed to eliminate excessive data 
computations and provide an accurate and efficient image reconstruction method. 

II. Relating Integrals of Line Integrals to a Distribution of Radioactivity 

In another embodiment, a plurality of integral of line integrals may be calculated on 
cone projections. For example, let the symbol Scb (*, ft \p) denote the integrated cone-beam 
line-integrals of the distribution of radioactivity on the one sheet cone whose apex is the point 
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axis of symmetry is the unit vector ft and a half angle is \f/. Assuming the vertex of the line 
integrals lies outside the distribution of the object, the cone-beam line integrals of the 
distribution may be defined as 



A 00 



g(<D,a)= J f(® + ta)dt. (21) 
o 

Using M T and a (</>,i// ) as previously defined in Equation (1) and Equation (3), the 
integrated cone-beam line integrals of the distribution can be defined as: 

A ^ n 

S CB (<D,^) = $ g((t>,M J cc(<t>,y,))d<{>. (22) 



The definition of the integrated cone-beam line-integral given in the Equation (22) and 
the surface integral given in Equation (2) may be compared for differences. First, there is a 
V radial weighting in the "dr" integral in Equation (2) that does not appear in Equation (22). 
This difference is significant because there is no simple way to obtain S ($, j8, \f/) from S CB 
($, ft or vice versa. Secondly, the surface integral given in Equation (2) there is a sin 
term that does not appear in the Equation (22). This leads to different values for these two 
integrals at yj, = 0. The surface integral at \J, = 0 equals zero. In this case, one data point may 
be known ahead of time. Furthermore, for the surface integrals, there may be an infinite 
number of significant digits. In contrast, at \J/ = 0 the integrated cone-beam line-integral may 
not always equal to zero. In fact, if the axis of symmetry associate with the integral intersects 
the distribution, then the integral is equal to 2tt times the line integral of the distribution along 
the axis of symmetry. 

Proof 3 

K 

F{p, <D • p) = lim f S CB (O, J3, y/)H s (cos y,) sin y, dy. (23) 

£ -° o 

Proof: 

Letting g(t) = H £ (d> p - 1) in Equation (7) and making changes to the variables results in 
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00 

J f{p,t)H e (0.fi-t)dt= $ /(<D - X )H e ( Z ■ P) d X . (24) 

-co 9? 3 

Letting ^ = ra for a e S 2 the right hand side of the above equation becomes: 

00 

JJ" f(®-ra)H £ (ra-/3)dar 2 dr. (25) 

OS 2 

Evaluating the above spherical integral in the same local coordinate system that was used in 
Equation (9) the above expression becomes: 

co2;r/r 

j J J / - rM T a(<f> y y/))H £ (r cos y/) sin <ty <ty r 2 dr. (26) 

0 0 0 

Now letting z = cos ^ results in 

oo 2^+1 

{ J J / (<D - rM T a(^cos-' z))H e (rz)dz d<j> r 2 dr. (27) 

0 0-1 

First, taking the limit as f->0 of left hand side of Equation (4) and the above expression, 
using Lemma 2 (discussed below) , and then letting y/ = cos -1 z in the resulting expression 
and exchanging integrals, and then using Equation (14) results in 

X 2 /TOO 

F(/3,® ■ 0)=lim j J J / (o - rM 7 a{<t>,v/))dr d<j> # £ (cos y) sin ^ dy. (28) 

e ~*° 0 0-1 

The conclusion of Proof 3 illustrates using Equations (21) and (22) in the right hand side of 
Equation (28). 



Lemma 2 

For a function g and a constant k 



Km \g(t)HXkt)dt = )im±- \ g{t)H £ (t)dt. 
Explicitly writing out the left hand side (LHS) of Equation (29) yields: 



= lim 



-1*1 



Replacing e/|^| with e results in 
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(29) 



(30) 



LHS = lim\ 

*-o k 



\-e b \-« e J 1 J 



(31) 



\-£ V-oo e J 

The conclusion of the Lemma 2 follows by observing that the right hand side of Equation (31) 
is equal to the right hand side of Equation (29). 

As such, a sequence of equations that relate the integrated cone-beam line-integrals to 
a distribution of radioactivity may be proven if the function F{pj) is known on a set P F , then 
the function/^) for all X e 91 3 may be found. The function F{fi,l) may be obtained on the 
set P F from values of S CB (<£,y9,^) via Proof 3 if the following condition is true: 

If for almost every (pj) e P F there exists an <D such that ®-j3 = £ where 
S CB (®,y!?,^) is known for all 0 < y/ < n . 

For convenience let the symbol P SCB represent any set of (<t>,j3,y) that satisfies the above 
condition. Using this symbol, this method of reconstructing the distribution of radioactivity 
from such a set of integrated cone-beam line-integrals can be representative symbolically as: 

S CB (0,/?,^)on Pscb F(PJ)on P F > /(^on 5R\ (32) 

As such, set Pscb and P F define what integrals are needed to eliminate excessive data 
computations and provide an accurate and efficient image reconstruction method. 

EXAMPLES 

Specific embodiments of the invention will now be further described by the following, 
nonlimiting examples which will serve to illustrate in some detail various features. The 
following examples are included to facilitate an understanding of ways in which the invention 
may be practiced. It should be appreciated that the examples which follow represent 
embodiments discovered to function well in the practice of the invention, and thus can be 
considered to constitute preferred modes for the practice of the invention. However, it should 
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be appreciated that many changes can be made in the exemplary embodiments which are 
disclosed while still obtaining like or similar result without departing from the spirit and scope 
of the invention. Accordingly, the examples should not be construed as limiting the scope of 
the invention. 

Example 1 

As mentioned above, the surface integral model and the integrals of line integral model 
each included a set of needed integrals that when reconstructed, may obtain an image. A 
completeness condition as defined herein, allows for determining what data is needed to 
complete this set. The surface integral completeness condition is: 

If on almost every plane that intersects a sphere with a radius bigger than the 
distribution, there is an apex where all surface integrals emanating from the apex 
whose axis of symmetry is normal to the plane are known, then the distribution of 
radioactivity can be obtained from the known integrals. 

Similarly, the integrated line integral completeness condition is: 

If on almost every plane that intersects the distribution, there is an apex where all the 
integrated cone-beam line-integrals emanating from the apex whose axis of symmetry 
is normal to the plane are known, then the distribution of radioactivity can be obtained 
from the known integrated line-integrals. 

In one embodiment, a cross-section of a sphere 200 contains a distribution of 
radioactivity as illustrated in FIG. 3. Sphere 200 has a radius 201 with a value of R + e. FIG. 
3 also shows a sine-on-the-cylinder trajectory 202. The sine-on-the-cylinder curve may be 
defined as a geometry that may include two periods of a sinusoid that may have been wrapped 
around a cylindrical surface. In particular, the geometry may be described by a vector-valued 
function O(A) = (c, cosA,c, sin X,c 2 sin 2,1) for 0 <X <2tt, where c, and c 2 may be selected 



25389513.1 



15 



so that a completeness condition is satisfied. (Smith, 1990). Suppose a first detector moves 
about the sphere 200 along trajectory 202. An arbitrary may intersect the trajectory 2, 3, or, 4 
times. For example, plane 203 intersects the trajectory 203 four times. A line {e.g., line 204), 
normal to plane 203 at the intersection of the sine-on-the-cylinder trajectory 202, must have a 
second detector that intersects the line 204 in order to satisfy the surface integral completeness 
condition. 

In an alternative embodiment, a detector may be large enough to subscribe a circle 
containing the distribution of radioactivity, where the circle has a radius larger than R + e. If 
the detector rotates along a circular trajectory about the sphere, then the surface integral 
completeness condition can be satisfied. In one embodiment, if the distribution is too large 
such that the detector can not subscribe the sphere, the detector may move along a sine-on- 
the-cylinder trajectory to satisfy the surface integral completeness condition. 

The completeness condition may further reduce the number of measurements taken. 
For example, for a given pair of detector elements, if a plane that intersects the first element 
and is perpendicular to the line that connects the two elements fails to intersect a sphere with 
radius of R + E , then the data associated with these two elements do not have to be processed, 
and hence do not have to be measured. 



The configuration of the detectors may also be considered to satisfy the completeness 
condition. In one embodiment, the surface integral completeness condition may be satisfied 
in which a first detector may be a planar detector and a second detector may be a spherical 
(e.g., hemisphere shape) detector. The detectors should be large enough such that a circle 
with a radius R + e may be subscribed on the detectors and a camera rotated along a circular 
trajectory about the ball that contains the distribution. 

In yet another embodiment, the first and second detectors may be spherical-shaped 
detectors. Again, if the detectors are large enough to embody a circle with radius R + e with 
a camera including the detectors rotate along a circular trajectory about the ball that contains 
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the distribution, the completeness condition is met. 

Alternatively, the first detector may be a cylindrical detector and the second detector 
may be a spherical-shaped detector. Furthermore, the first and second detectors may be planar 
detectors. Further, there may be more that two detectors used to obtain conical data to satisfy 
the completeness condition. 



Example 2 

The Compton data is considered to be integrals, such as surface integrals or integrals 
of line integrals, over the points that lie on a cone that have been weighted by the Klein- 
Nishina distribution of scatter angles and blurred by the angular Doppler broadening. 
Furthermore, the data may be significantly random because of inherent randomness of the 
generation of photons. To develop a computational efficient reconstruction method that will 
mitigate the effects of this weighting, blurring, and randomness, a two-step image 
reconstruction method may be implemented. The first step may mitigate the effects of the 
weighting, blurring, and randomness to obtain a good estimate of conical surface integrals. 
The second step may perform a computational efficient tomographic reconstruction from 
estimates of the conical integrals obtained in the first step. 

Step 1 

First, the energy that is measured may be partitioned into "energy bins" For k = 1, - , 
N e , the parameters e k may be chosen such that e k < e k+ i and 0 ^e k <e max for all k, where emax 
may be the maximum energy level to be measured. The interval [e k , e k+ i] forms the k th energy 
bin. Next let Y (j, 1, k) be the random variable associated with the number of photons that 
interact with the j th first detector and the 1 th detector and are counted in the k th energy bin. 

As mentioned previously, the purpose of the first step of the reconstruction method is 
to mitigate the effects of the weighting, blurring, and the randomness. This mitigation may be 
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performed independently for each pair of detectors in the first and second detector. For a 
given first detector j' and a given second detector 1' the following may be defined. First the 
reconstruction volume is partitioned into hollow cone regions. As shown in Figure 1, the 
variable yp is the angle measured from the line that connects the intersection point ($ and 
of a photon onto the two detectors. For m = 1, - , N the ^ m 's are chosen such that \p m < ^ m+ , 
and 0 ^ m &v for all m. The selection of the ^ m ' s are not necessarily depended upon the 
selection of the eic's. For example, geometrically each ^ m defines a cone whose apex is at the 
center of the first detector and axis of symmetry is the line that connects the centers of the 
detectors and the "half width" is the angle The m th hollow cone region is defined to be 
the portion of the reconstruction volume that lies in between the \p m and ^ m+1 cones. It should 
be noted for future reference that this partition is dependent upon j' and 1'. 

Having defined the term "hollow cone region" it can be explained why angular 
Doppler broadening may be thought of as having the effect of blurring the data. In the 
absences of angular Doppler broadening all the photons that loss a certain energy level would 
originate in one hollow cone region. This region would be determined by Compton's law. In 
the presents of angular Doppler broadening these photons, generally speaking, would originate 
in more than one region. This spreading out of the photons (i.e., the blurring of the photons) 
will degrade the quality of the images that result. Unfortunately, this blurring is shift variant 
which makes mitigating its negative effects difficult to do. [Evans et ah, 1999] 

Let W(j', P, k, m) be the random variable associated with the photons generated in the 
m th hollow cone region that are counted in the k th energy bin and interact with the fixed 
detector pair. Thus it can be written: 

Y(j\P,k) = £«OV\*,*)fork=l;-;N. (33) 

m=l 

Assume that the number of photons generated in the regions are independent and 
are Poisson distributed with mean o(j', P, m). The probability p(j', 1', k | m) is defined to be 
the conditional probability that a photon will interact with the fixed detector pair and will be 
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counted in the k th energy bin given that it was generated in the m th hollow region. From 
Equation (1) and from the properties of expected values of Poisson random variables [Feller, 
1968] it can be written; 

E {Y (j',1', k)} = XpU'J'J I m)a(fl',m) for k = 1; •■ ;N e (34) 

Note from Equation (33) that Y(j\ 1', k) is a sum of independent Poisson random 
variables and hence is a Poisson random variable. As such, the maximum likelihood estimate 
of Y(j\ P, k) is the number of photons measured in the k th energy bin for the fixed detector 
pair. This number will be denoted as n(j\ P, k). Substituting the estimate in for the right 
hand side of Equation (34) yields: 



n(j', P, k) = Zp(J',l',k | m)a{f,V,m) ; for k = 1; - ;Ne 



(35) 



By using an algorithm that is known in the Art, such as a Penalized Weighted Least 
Square or a ML-EM algorithm, a Q\ P, m) for m = 1; ••• ;N may be determined. The values of 
ot(j\ P, m) will be used in Step 2. 

Equation (35) has a number of very advantageous features from the computational 
point of view. First, the processing for estimating the a (j\ P, m)'s may be done 
independently for each pair of detectors which allows for the reconstruction method to be 
parallelizable. Namely, all the processing for a given detector pair may done on a given node 
of a parallel computer. This may allow for a reduce execution time if a parallel processing 
computer is used. Secondly, the values of p(j', P, k | m)'s for k = 1; -;N e and m = 1; - ;N 
may be stored on a given node. If a 51 1 keV isotope is used and if the width of the energy bin 
is about 1 keV then N e «10 +2 - If Ni/> «10 +2 as well, then the storage requirement per node 
would be «10 4 , which is an acceptable number. Hence the probability values can be pre- 
computed, storage, and read into the program when needed. Furthermore, each node would 
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have the probability values it needs and no probability value would have to be shared from 
node to node and thus improving the effectiveness of the algorithm. 



To find an appropriate expression for calculating p(j', 1', k | m) for Equation (34), it 
may be assumed that the amount of energy loss by a photon is independent of the fact that the 
photon interacted with the j' th first detector and is also independent of the fact that it 
originated in the m th region. Using the definition of conditional probability, it can then be 
written that 

p(j', 1', k | m) = p(l' | j', k, m) pG' | m) p(k) (36) 

where 

■ p(l' | j', k, m) is the probability that a photon will interact with the l' th second detector 
given that it interacted with the j' th first detector and was counted in the k th energy bin and 
originated in the m th region; 

■ p(j ' | m) is the probability that a photon will interact with the j' th detector given that the 
photon was generated in the m th region; and 

■ p(k) is the probability that a photon will be counted in the k th energy bin. 

To find an appropriate expression for p(l\ j', k | m) Equation (36) needs to be 
modified. To obtain a relatively accurate value for p(l', j', k, m) the m th hollow cone region is 
partitioned into sub-regions. First the angle f may be defined as an angle around the line 
connecting the centers of the detectors. The domain of which is [0, 2x], may be divided 
into N t equal length intervals. Additionally, the interval W m , ^ is divided into Ns equal 
length intervals. The (s,t) th sub-region of the m th hollow cone region is the region defined by 
the interval [ft, and the interval [* m + (i - 1)( ^ m+] - ^ m)/Ns , ^ + (i)( ^ . 

Using the definition of conditional probability and the total probability theorem, it can be 
written that 
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i;i>(r|/,*,*,OKAiUO 

P(l' I J", k, m) = — — (37) 

p(f,k,m) 

where 

■ PO'I j', k, s, t) is the conditional probability that the photon will interact with the l' th 
second detector given that it interacted with the j' ,h first detector and was counted in the 
k th energy bin and originated in the (s, t) th subregion; 

" PG\ k, s, t) is the probability that a photon will originate in the (s, t) th sub-region and will 
be counted in the k th energy bin and will interact with the j' th first detector; and 

■ p(j\ k, m) is the probability that a photon will originate in the m th region and will be 
counted in the k th energy bin and will interact with the j' th first detector. 

Since the sub-regions are assumed to be small and are all within the same region, the 
PC, k, s, t)'s are assumed constant with respect to the s and t. Assuming as before that the 
event k is independent of the other events, then the above equation can be written as 



Now Equation (36) can be written as: 

PG \ 1\ k |m) = p(k) p(j' | s, t) ^4£f>(/l fXs >t ) (39) 

P\ m ) <=i s=i 

where 

■ p(j ' | s, t) is the conditional probability that a photon will interact with the j ' th first detector 
given that the photon originated in the (s, t) th sub-region; 

■ p(s, t) is the probability that a photon was generated in the (s,t) th sub-region; and 

■ p(m) is the probability that a photon was generated in the m th region. 

To find an appropriate expression for the p(k), the "Klein-Nishina distribution," m 
&»(•), may be integrated, namely: 

21 

25389513.1 



p(k)= t' h(g)f KN (g)dg (40) 

where the function h(g) accounts for any variation in the detector efficiency. The probability 
pG'l s, t) may be calculated by performing a numerical integration of the steroangle subtended 
by the j' th detector from each point inside the (s, t) th sub-region. To compute the ratio of the 
probabilities p(s, t)/p(m) it is assumed the point within the m ,h region at which a photon was 
generated is equally likely. Because of this, the ratio will be taken to be the volume of the 
(s,t) th sub-region divided by the volume of the m th region. 

An appropriate expression for p(l' | j', k, s, t) for Equation (39) may found by 
integrating a joint probability density function. The angle 6, as illustrated in FIG. 4 which is 
the scatter angle, is measured from the center ray of the (s,t) th sub-region. The angle 0, which 
is the angle "around" this line, is not shown in FIG. 4 for the sake of berevity. 

To define a joint probability density, it is assumed that the <f> angle associated with the 
scatter of the photon is independent of the angle d. This allows the joint probability density 
function to be written as a product of marginal probability density functions. It is further 
assumed that <\> is uniformly distributed on the interval [0, 2tt]. The marginal probability 
density function associated with d is taken to be the function f D (., e ), the angular Doppler 
broadening at the fixed energy level. Thus, the joint probability density function f(0, $) for the 
k th energy bin may be defined as 



f(d, <t>) = j h™^?*- ^x0<e<nand0<<p<2n; 

|0 otherwise ( ' 



An appropriate expression for p(l' | j', k, m) may be found by integrating the above 
joint probability density function. By defining « = (j\ 1', s, t) to be the set of ordered pairs 
of (6, eft that result in the photon being detected in the l' th detector, it can be written that; 



25389513.1 



22 



p(i'|j\k, M )= ll^^jv^deicp 



(42) 



In another embodiment, an alternative means to using Equation (39) may be values 
obtained empirically for the needed probabilities. For example, by placing a small radioactive 
source inside various locations with the m th region, the p(j', 1', k | m)'s for the various values 
of k may be taken to be a relative frequency of an observation. In yet another alternative, 
these relative frequencies may be obtained via Monte Carlo computer simulations. An 
advantage of this latter method is that self absorption may be taken into account. 

Step 2. 

The a (j\ 1', m)'s obtained in Step 1 may be an input to the reconstruction method 
used in this step. In one embodiment, surface integrals models or the integrals of line integral 
models, as described above, may be used in the reconstruction method. For example, 
Equation (20) and (32) may be used to reconstruct an image. Alternatively, an ART-like or a 
SIRT-like algorithm may be developed. Further, an ML-EM may be implemented. For 
example, by voxelizing the reconstruction volume, an equation analogous to Equation (35) 
may be developed. However, the probability term would be a conditional probability that a 
photon will interact with the fixed element pair and will be counted in the k th energy bin that 
may be generated for a certain voxel. 

* * * 

With the benefit of the present disclosure, those having skill in the art will comprehend that 
techniques claimed herein and described above may be modified and applied to a number of 
additional, different applications, achieving the same or a similar result. The claims cover all 
modifications that fall within the scope and spirit of this disclosure. 
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